Velocity distribution in granular gases of viscoelastic particles 
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The velocity distribution in a homogeneously cooling granular gas has been studied in the vis- 
coelastic regime when the restitution coefficient of colliding particles depends on the impact velocity. 
We show that for viscoelastic particles the simple scaling hypothesis is violated, i.e., that the time 
dependence of the velocity distribution does not scale with the mean square velocity as in the case 
of particles interacting via a constant restitution coefficient. The deviation from the Maxwellian 
distribution does not depend on time monotonously. For the case of small dissipation we detected 
two regimes of evolution of the velocity distribution function: Starting from the initial Maxwellian 
distribution, the deviation first increases with time on a collision time-scale saturating at some max- 
imal value; then it decays to zero on much larger time-scale which corresponds to the temperature 
relaxation. For larger values of the dissipation parameter there appears an additional intermediate 
relaxation regime. Analytical calculations for small dissipation agrees well with the results of a 
numerical analysis. 



PACS numbers: 81.05.Rm, 36.40. Sx, 51.20.+d, 66.30.Hs 



I. INTRODUCTION 

The statistical properties of granular gases have been 
intensively studied in recent time, in particular with re- 
spect to cluster formation process, e.g. [jjj and other- 
structure formation, e.g. In the present paper we 
are concerned with dynamical processes in granular gases 
which precede clustering, i.e. in the homogeneously cool- 
ing state (HCS). In difference to the state when particles 
form clusters and other long range structures, in the HCS 
(due to its definition) one may drop the explicit spatial 
dependence of the statistical properties, which simplifies 
an application of standard methods of the gas kinetic 
gas theory. Granular gases in the HCS were intensively 
investigated recently (see e.g. || for a review) focusing 
on the velocity distribution function which is one of the 
most important characteristics of the system of granular 
particles. It has been argued that the distribution func- 
tion might deviate from the Maxwellian 0,0, and this 
deviation has been also quantified dUQ] • 

In all of these studies a constant restitution coefficient, 
characterizing the energy loss due to a particle collision 
was assumed. The restitution coefficient relates the ve- 
locities of the colliding particles before a collision V\ , V2 
to the velocities after the collision v{ , : 

1 



v{ = ui - -(1 + e)(ui2 • e)e 
v% = v 2 + -(1 + e)(v 12 ■ e)e 



(1) 



where vi 2 — v\ — v 2 is the relative velocity and the unit 
vector e — r 12 / Iri^l gives the direction of the inter-center 
vector fi2 — r\—T2 at the instant of the collision. Strictly 
speaking the restitution coefficient e as introduced in 



Eq. (Q) describes the collision of smooth inelastic par- 
ticles, when only the normal component (vi 2 ■ e) of the 
relative velocity v\ 2 changes. Therefore, it is termed as 
normal restitution coefficient. Using the tangential resti- 
tution coefficient |^,|8|,[l0| , one can account for the change 
of tangential component of the relative velocity at the 
collision of rough inelastic particles. In what follows we 
assume that the particles are smooth and the dynamics 
of a collision is completely described by the change of the 
normal component of the relative velocity. 

Experiments, as well as theoretical studies show, how- 
ever, that e noticeably depends on the impact velocity 
v\2 Jl2]-|l4|]; even a dimension analysis shows that the 
assumption of the constant restitution coefficient contra- 
dicts physical reality Jl5],[l6| . This dependence may cause 
rather important consequences for various problems in 
granular gas dynamics Jl7],[l8| • The problem of the resti- 
tution coefficient's dependence on the impact velocity has 
been addressed in [p||ll|] , where the generalization of the 
Hertz contact problem was developed for the collision 
of viscoelastic particles (scaling analysis of this depen- 
dence has been also addressed in fl3|j) . The generalized 
Hertz collision equation derived in H has been solved 
analytically to obtain the velocity-dependent restitution 
coefficient |l^] 

e = 1 - dAa 2 / 5 \v 12 ■ e1 1/5 + C 2 A 2 a i ' b \v 12 ■ ef /5 T ■ 

(2) 



with 



3/2 



TO eff (l - V 2 ) 



(3) 



where Y is the Young modulus, v is the Poisson ratio, 
R cS = RiR 2 /{Ri + R 2 ), m cS = rmm 2 /(mi + m 2 ) (R 1/2 
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and mi/2 are radii and masses of colliding particles), A is 
the dissipative constant, which depends on the material 
parameters (see |8| for details). Numerical values for the 
constants C\ and Ci obtained in jl9| may be also written 
in a more convenient form |16|: 



d = 



C'2 



r(3/5)V? 

2 l/5 5 2/5 r (21/l ) 



= 1.15344, 



(4) 
(5) 



Although the next-order coefficients C3 = 0.315119CJ 1 , 
C 4 = 0.161170C1 1 , are now available @, we assume that 
the dissipative constant A is small enough to ignore these 
high-order terms. 

The aim of the present study is to analyze how the 
impact-velocity dependent restitution coefficient given by 
Eq. (||) for the collision of viscoelastic spheres, influences 
the velocity distribution in a granular gas of identical par- 
ticles in HCS. To address this problem we use the Sonine 
polynomials expansion for the velocity distribution func- 
tion and analyze the time-dependence of the expansion 
coefficients. 

In Sec. we introduce the necessary variables, briefly 
sketch the method of Sonine polynomial expansion and 
summarize the knowledge about the velocity distribution 
function in granular gases under the assumption of a 
constant restitution coefficient. In Sec. [II we analyze 
the Boltzmann equation for the granular gas with the 
velocity- dependent e in the HCS and calculate the first 
few coefficients of the Sonine polynomials expansion. We 
show that these coefficients occur to be time-dependent, 
so that the velocity distribution function does not have 
a simple scaling form. In Sec. ^ we consider the time 
evolution of temperature and of the velocity distribution. 
The high- velocity tail of the distribution function is an- 
alyzed in Sec. |V[ In Conclusion we summarize our find- 
ings. Some technical detail of the calculations are given 
in the Appendixes. 



II. SONINE POLYNOMIAL EXPANSION FOR 
GRANULAR GASES 



m is the mass of the granular particles, and d is the di- 
mension. The temperature is related to the second mo- 
ment of the velocity distribution in the same way as for 
equilibrium molecular systems: 



mv 

dv——f{v,t) 



(8) 



Then the expansion of the scaling function f{c) (where 
c = v/vo(t)) in terms of the Sonine polynomials reads 



(9) 



where 4>(c) = n~ d / 2 exp(— c 2 ) is the Maxwellian distribu- 
tion for the rescaled velocity. The Sonine polynomials 
S p (c 2 ) satisfy the orthogonality conditions 



dc(j>(c)S p (c 2 )S p >(c 2 ) = 5pp>Afp 



(10) 



with 5 P p> being the Kronecker delta and with the normal- 



ization constant Af p j4]||]. For dimension d = 3 which is 
addressed here the first few Sonine polynomials read 

S (x) = 1 

Sl (x) = -x 2 + l (11) 



S 2 (x) 



-x 

,•2 



bx 15 
T + ~8~ 



(12) 



The coefficients a p of the expansion may be found as the 
polynomial moments of the function /(c) 



(13) 



The coefficients a p do not depend on time for a constant 
restitution coefficient [E0|. These were first applied for 
the granular gas in Ref. |4] and then recalculated recently 



For granular gases where the particles interact via a 
restitution coefficient e = const it was argued that the ve- 
locity distribution f{v,t) has the scaling form, i.e. that 
its time-dependence may be written as (here we follow 
notations of Ref. H ) 



Mt) 



-f 



(6) 



where n is the number density of the granular gas, vo(t) is 
the thermal velocity, defined in terms of the temperature 
of the granular gas 



T(t) = \mv 2 {t) 



(7) 



ax = 



a 2 



16(l-e)(l-2e 2 ) 



9 + 24d + 8ed + Ale + 30(1 - e)e 2 



(14) 
(15) 



The first relation (|14|) follows from the definition of the 
temperature of the granular gas (this we explain in more 
detail below), while Eq. (15) has been obtained within 
the linear approximation with respect to 02 . Complete 
analysis, which goes beyond the linear approximation, 
has been performed [JtJ , and it has been shown JjJ that 
the linear solution ( |15[ ) is rather accurate for the whole 
range of e with a maximal deviation from a total one 
less than 10% [^lj. All the higher-order coefficients were 
neglected under assumption of small deviations from the 



2 



Maxwellian distribution. Since a p do not depend on time, 
the scaling form of the velocity distribution function (|^) 
persists with time for the case of e = const. 

Since the average velocity of a granular gas decreases 
due to permanently decreasing temperature, the "typi- 
cal" restitution coefficient will increase with time as it 
follows from Eq. (||). Thus one can expect that the co- 
efficients of the Sonine polynomials expansion, which de- 
pend on the restitution coefficient (see e.g. (|l5j)) should 
also change with time. This conclusion, however, contra- 
dicts the assumption, that the scaling function (||) docs 
not depend on time and implies that the common scheme 
of calculation of the Sonine polynomials expansion coeffi- 
cients breaks down if e is not a constant. For for the latter 
case, one needs to develop a more general approach. 



III. KINETIC EQUATION FOR THE 
COEFFICIENTS OF THE SONINE 
POLYNOMIALS EXPANSION 

We start from the Enskog-Boltzmann equation for the 
distribution function /(r, v, t) for a granular gas of inelas- 
tic spheres which in the force-free case does not depend 
on r, Hence, one can write flp2| 



92(a)! (f J) 



(16) 



d_ 



f(v\,t) = g2(o)o 2 J dv 2 J de&(-v 12 ■ e) |« 12 • e\ x 



where a is the diameter of the particles. The contact 
value of the two-particle correlation function, 52(c) = 



(2 — 7j)/2(l — rj) 3 ]23| (with V = \ ^na 3 being the packing 
fraction) accounts for the increasing collision frequency 
due to the excluded volume effects. 9(x) is the Heavi- 
side step-function. The velocities v*{* and v^* refer to the 
pre-collisional velocities of the so-called inverse collision, 
which results with v\ and v 2 as the after-collisional ve- 
locities (the relation between these velocities are similar 
to that of Eq. (Q), but with the impact- velocity depen- 
dent restitution coefficient, see Appendix A). Finally the 
factor 



1 + ^dAa 2 ^ 5 |0i2 • el 
5 

g^AV/ 5 \v 12 ■ ef/ 5 



1/5 



(17) 



in the gain term appears respectively from the Jacobian 
of the transformation dv[*dv^* — > dv\dv 2 and from the 
relation between the lengths of the collisional cylinders 
£ \vi2 ■ e\dt = \vi 2 ■ e| dt (see Appendix A for details). For 
the constant restitution coefficient x = l/ e2 = const 

Some important properties of the collisional integral 
hold also for the case of the impact-velocity dependent 
restitution coefficient. Namely, it may be shown that the 
relation 



d f c) 

j t m)) = / dv 1 i>(v l )—f(v u t)= I dvnP(vi)I(f,f) 



(18) 



g 2 {a)cr 2 



dvidv 2 / de®(-vi2 ■ e) |€T 12 • e] f(vi,t)f(v 2 ,t) A [ip(vi) + ip(v 2 )] 



holds true, where (ip(t)) = J dviip(v)f(v,t) is the aver- 
age of some function ip(v), and Aijj(vi) = [ip(v* ) — ip{vi)] 
denotes change of ip(i>i) in a direct collision. 

Now we analyze the scaling ansatz (^) for the velocity 
distribution function. Using this ansatz and performing 
calculations similar to that in Ref. ||, one would find 
corresponding expressions for the coefficients a p of the 
Sonine polynomial expansion. These would occur to be 
time-dependent due to permanently decreasing average 
velocity of the cooling gas and thus permanently increas- 
ing effective value of the restitution coefficient. This how- 
ever means that the simple scaling (^) for the velocity dis- 
tribution function does not hold for the case of interest. 
Technically, as we show below, this follows from the ad- 
ditional timc-dcpcndcncc of the factor \ m the collisional 
integral, which does not depend on time for e = const. 

Thus, it seems natural to write the three dimensional 
distribution function in the general form 



f(v,t) 



f(c,t) 



(19) 



f(c,t) = cf>(c)ll + J2a P (t)S p (c 
{ P =i 



(20) 



with 



and find then equations for the time- dependent coeffi- 
cients a p (t). 

Substituting (|l^) into the Boltzmann equation ( |l6| ) we 
obtain 

1 dvr, / d \ ~,_ . 1 d ~,_ . 

~~TT [ 3 + Cl 7T~ /( ? i>*) + — ~7/(ci,i) = 
Vq dt \ da J vo ot 

= .g 2 (a)a 2 n/(/,/) , (21) 
where we define the dimensionless collisional integral 

I (/, /) = J dc 2 J deQ(-c 12 ■ e) |c 12 • e\ X 

x {x/(cT,i)/(?,t)-M,t)/(4i)) • (22) 
The reduced factor \ 



3 



1 i ^nx>\- "i 1 / 5 i 66 2 /2 2/5 
X = 1 + Y 1 ' 12 ' ^ 25 1 |ci2 ■ e| 



(23) 



depends now on time via a quantity 

6'(t) = Aa 2 ' 5 [2T{t)] 1/w = 6 [2T(i)/T ] 1/10 . (24) 

Here S = Aa 2 ^ 5 (Xb) 1 ^ 10 , To is the initial temperature, 
and for simplicity we assume the particles to be of unit 
mass, m = 1. 

The rate of change of the thermal velocity dvo / dt in 
( pl| ) may be expressed in terms of the temperature decay 
rate dT /dt, which reads according to the definitions 
and relation (|l8|) for the time derivatives of averages 



dT 1 



:9i 



{a)a 2 nvl J dc lC 2 I (/ r , /) = -^BT m . (25) 



We define here B = B(t) = v (t)g 2 (a)a 2 n and introduce 
the moments of the dimensionless collision integral: 

th = ~ J d& 1( fj(fj) . (26) 
With this notations we recast (Ell) into the form: 



(27) 



f (3 +C1 a) m +B -^ t m =/(/./ 



Note that contrary to the case of e = const, where x = 
1/e 2 = const, the factor x depends now on time, which 
does not allow to write the collision integral in terms 
of only scaling variables. This implies time dependence 
of all the moments [i v (which were time-independent 
for constant restitution coefficient) and correspondingly 
causes time dependence of the Sonine polynonomials ex- 
pansion coefficients a p . 

Multiplying both sides of Eq. (|27]) with c\ and inte- 
grating over ci we obtain: 



^-p (c p ) - B 1 ^2 h kV kp = A* P 1 
6 fc=i 



(28) 



where integration by parts has been performed and where 
we define 



Vk P = / 4>{c)c p S k (c 2 )dc 



and 



c p f(c, t)dc. 



(29) 



(30) 



The calculation of the Vkp is straightforward; the first few 
of them read: 



n 15 

v 2 2 =0; 2/ 24 = — . 



(31) 



The odd moments /c 2rl+1 ) are zero, while the even ones, 
(c 2n '), may be expressed in terms of with < k < n. 
This follows from the fact that c 2n may be written as a 
sum of Sonine polynomials Sk{c 2 ) with < k < n and 
from the orthogonality condition (JTo|) . Namely, using 
c 2 = |S'o(c 2 ) — Si(c 2 ) together with the expansion ( |20| ) 
and condition (|l0|), one easily finds 



<c 2 > 



dc<j)(c) 



ai 

2 2 



-S* (c 2 ) - S^c 2 ) 



,fc=o 



(32) 



with ao = 1 and where we use the normalization constant 
Wi = I [see Eq. (Jl~0|) ] . From the definitions of temper- 
ature and of the thermal velocity (|7), (||) follows that 
(c 2 ) = § (see also ©). Then Eq. (§!) implies a x = 
in accordance with Rcf. B . Similar considerations yield 



15 



(I + 02) 



(33) 



The moments /i p may be also expressed in terms of coef- 
ficients CI2, 0,3, ■ ■ •, therefore, the system (28) is an infinite 



(but closed) set of equations for these coefficients. 

It is not possible to get a general solution to the prob- 
lem. However, since the dissipative parameter S is sup- 
posed to be small, the deviations from the Maxwellian 
distribution are not presumably large. Thus we assume, 
that one can neglect all the hi gh-o rder terms in the ex- 
pansion ( pcf ) with p > 2. Then (p8|) is an equation for the 
coefficient a 2 . For p = 2 Eq. (|28j ) converts into identity 
since (c 2 ) = |, ai = and due to (fji"|). For p = 4 we 
obtain 



4 4 
d 2 - - Bp 2 (1 + a 2 ) + —-8/^4 = . 
3 15 



(34) 



where the relations (|31j) and ( |33| ) have been used. In Eq. 
B depends on time as 



B(<) = (8^)- 1 / 2 t c (0)- 1 [T(<)/T ] 1 / 2 



(35) 



where To is the initial temperature and r c (0) is the initial 
mean-collision time 



M- 1 =^ 1/2 92{<r)a 2 nT 1 J 2 



(36) 



The time evolution of the temperature is determined by 
Eq. (|25|), i.e., by the time dependence of /.i 2 . 

The time-dependent coefficients p p (t) may be ex- 
pressed in terms of a 2 according to definition ( p6| ) and 
the approximation / = </>(c) [l + a 2 (t)S 2 (c 2 )] . One ob- 
tains: 

M P = J J d ^2 J deO(-ci2 ■ e) |c 12 • e| (f>(cx)<p(c 2 ) 

x{l + a 2 [S 2 (c 2 ) + S 2 (c 2 )] + a 2 S 2 (c 2 )S 2 (c 2 )} x 
xA( c p + c^) (37) 



4 



with the definition of A(c^ + c^) given above. After long 
and tedious calculations (details are given in Appendix 
B) one arrives at the following result for the moments: 

M2 = 5' [Ai + A 2 a 2 +A 3 al] - s ' 2 [&L + A 5 a 2 + Asa 2 ,] 

(38) 

and 

M4 = [Bi + B 2 a 2 + B 3 a 2 2 ] + 5 ' [B 4 + B 5 a 2 + B 6 a 2 2 ] (39) 

-5' 2 [B, + B 8 a 2 + B 9 a 2 ] 

where A n and B n are pure numbers. The coefficients A n 
read 



A\ — uJq A 2 — £rOJo A3 — oifin^O 



25 
119 



2500 L 
4641 



A4, - L>! A 5 = 400^1 A 6 - 7340000 Wl 



(40) 



with 



(41) 



uj = 2V2^2 1/10 r Ci = 6.48562 . . . 

w x = V2^2 1/5 T fy J C\ = 9.28569. .., (42) 



and the coefficients B n are 



S4 = yq 



B 2 
B 5 



4V27T 
1806 



250 



w 



£ 7 = Swi B 8 = 



149054 
13750 



#3 = 
#6 = 
#9 = 



1^ 
567 



12500 



0JQ 



(43) 



348424 
5500000 



Thus, Eqs. (B4j) and (|25j), together with Eqs. (|24|), 
(|35|), ( |3"s| ) and p9j ) form a closed set to find the time 
evolution of the temperature and coefficient a 2 . We want 
to stress an important difference for the time evolution 
of temperature for the case of the impact- velocity depen- 
dent restitution coefficient, compared to that of the con- 
stant restitution coefficient. In the former case it is cou- 
pled to the time evolution of the coefficient a 2 , while in 
the latter case there is no such coupling since a 2 — const. 
This coupling may lead in some case to rather peculiar 
time-dependence of the temperature. The problem of the 
time dependence of temperature and the velocity distri- 
bution function will be discussed in detail in the following 
section. 



IV. TIME EVOLUTION OF TEMPERATURE 
AND OF THE VELOCITY DISTRIBUTION 
FUNCTION 

To analyze the time evolution of the temperature and 
of the coefficient a 2l characterizing the velocity distri- 
bution function, we introduce the reduced temperature 
u(t) = T(t)/T and recast the set (§3), (§f) into the form 



u + r n u 



-.a 2 



500 



-1 r 17/10 119 1547 

'oV* 17 ' 10 (3 + 2lo a2 + i280oo< 

a 2 - tqu 1 / 2 ^ (1 + a 2 ) + -r u^ 2 ^4 = . 

5 



The characteristic time 
_i 16 



qv& ■ Tc(O)- 1 



(44) 
(45) 

(46) 



describes the time evolution of the temperature (see be- 
low), with 

9o = 2 1 / 5 r(21/10)Ci/8 = 5- 2/5 % /r?r(3/5)/8 

= 0.173318... (47) 
2 

TO 



T-c(O) 



-1 



3V2tt 

qt = 2 1 / w (lu 1 /uj q ) = 1.53445. 



(48) 
(49) 



To obtain these equations we use the expressions for 
/J. 2 (t), B(t) and for coefficients A n . Note that the char- 
acteristic time To is <5 _1 ^> 1 times larger than the mean 
collision time t- c (0). 

We will find the solution to these equations as expan- 
sions in terms of the small dissipative parameter 5 (see 
Eq.(H): 



u = u + S ■ Ul + S 2 ■ u 2 + ■ ■ ■ 
a 2 = a 2a + S ■ a 21 + 5 2 ■ a 22 + ■ 



(50) 
(51) 



Substituting Eqs. (f30|), (|5lJ), Q38p and fl39J) into Eqs. 
(0), (@), one can solve these equations perturbatively, 
for each order of 5. Collecting terms of the order of 0(1) 
we obtain: 



wo + t 



— — aoo 
3 5 u 



7 

500' 



a 2o ) u o /5 = 



a 20 + rwl' 2 I a 20 + —020 J = 



1 



32 



where 



n 



:r n B 2 



15 



Tc(O)- 1 



(52) 
(53) 

(54) 



and we use the definition of ro, and expressions ( ftlf ) for 
B 2 and B 3 , which are zero-order coefficients in the expan- 
sion of /Z4 on 5. Changing variables 



t -> r 



dt'u /2 (t') 



(55) 



in Eq. (|53|) one finds the solution of this (Riccati) equa- 
tion: 



a 20 {t) 



a 20 (0) 



[1 + ^0(0)] e T - ^om(0) 



(56) 



5 



According to Eq. (E2J) the characteristic time scale for 
uo(t) is t ^> r c (0), therefore, for t ~ t c (0) <C t one can 
approximate = T(t)/To « 1. Moreover, if the initial 
deviation from the Maxwellian distribution is not large, 
i.e. a2o(0)/32 <C 1, one can approximate for this time 
interval: 



O20 (*) « a 2 o(0)e- 4t/5TE(0) 



(57) 



with te = |r c being the Enskog relaxation time. There- 
fore, a,2o(t) vanishes for t ~ tq 3> r c (0). This refers to the 
relaxation of an initially non-Maxwellian velocity distri- 
bution to the Maxwellian one. Note that the relaxation 
occurs within few collisions per particle, similarly to the 
relaxation of molecular gases. 

We now assume that the initial distribution is 
Maxwellian, i.e., that O2o(0) = for t = 0. Then the de- 
viation from the Maxwellian distribution originates from 
the inelasticity of the particle collisions. For the case 
fl2o(0) = (and thus aa{t) = 0, see Eq.([56|)) the solution 
to Eq. rt52) reads 



u {t) 



To 



t 



-5/3 



(58) 



which coincides with the time-dependence of the temper- 
ature obtained previously using scaling arguments ]l9| ] 
(up to a constant To which may not be determined by 
scaling arguments). 

For the order 0(S) we obtain: 



8 



3/5 



8/5 



Hi + — u ' ui + - — u a 2 i 

OTq OTo 



3r 



17/10 n fKM 

qiu Q = (59) 



with 



ri 



O21 + riu^ 2 a 2 i + r 2 u^ 5 = (60) 



(J^ 2 1 /io ( 8rr)- 1 /2 {Bi - 5^) T -\0). (61) 



For i C To we have uq ~ 1 and Eq. ( p0| ) reduces to 

021 + ri02i = —r% (62) 
with the solution: 

Oai(t) = ~ (1 - e- rit ) = -h (l - e -«/6™(0)) (63) 

where 

h = ra/ri = f ^) T H±j 2 1 / 5 C 1 = 0.415964 (64) 

and we used the definitions of fj , r 2 and the values of „4 fc 
and Bk given above. As it follows from Eq. (|63|), after 
a transient time of the order of few collisions per parti- 
cle, i.e. for r E (0) < t <C t , a 2 (i) saturates to the value 
a 2 = —hS = —0.415964(5, i.e. it changes only slowly on 
the time-scale ~ r c (0). 



For t > To the rescaled temperature varies uo 
(t/To)- 5/3 [see Eq. ©], and Eq. @ reads 



021 + ri (V T o) 5/6 021 = -^2 {t/ro) 
Using the power-law ansatz 

021(4) ~ (t/T )- U 



(65) 



(66) 



the asymptotic analysis of Eq. (|65J) yields the exponent 
v = 1/6 and an estimate for the prefactor. Thus, we find 
for t > r : 



a 2 i (t) 



(t/ro) 



-1/6 



-h (t/r ) 



-1/6 



(67) 



Therefore, a 2 ±(t) decays to zero on the time-scale ~ To, 
i.e., slowly on the time-scale ~ t c (0) -C to. The ve- 
locity distribution, thus, tends asymptotically to the 
Maxwellian distribution. 

One can also find the general solution of Eq. (pCf): 



021 (t) = -6T r 2 exp|-6T ri (1 + £/t ) 1/6 | 

— dx . 

6rpri ^ 



Noticing that 



6t ?"i = (q 5) 1 

fi 12 A-1 
DT r 2 = —6 



(68) 

(69) 
(70) 



due to the definitions or r\, r 2 and To, one can write for 
a 2 (t) = S ■ a 2 i(t) in linear with respect to S approxima- 
tion: 



oa(t) = -^wity 1 {Li [w(t)} - Li HO)]} 
5 



where 



w(t) = expUqoS)- 1 (l + t/T f /6 



(71) 



(72) 



and Li(x) is the logarithmic integral. It is not difficult to 
show that from the general expressio n (7j ) both limiting 
dependencies ( |6^ ) for t <C t and ([57 ) for t 3> t are 
reproduced. 

We could not find the general solution for ui(t), how- 
ever, one can obtain the solution for t To- Substi- 
tuting asymptotic expressions u (t) ~ (£/t )~ 5 / 3 and 
o 2 i (i) ~ -h(t/T y 1/6 into Eq. @ for m(f) we recast 
this equation into the form 



8 , , x-i /2, 5 
ui + ~(t/T ) mi = I -o + -gi 



(4/to) 



-17/6 



(73) 



Again a power-law ansatz Oi(t) ~ (t/T ) a allows to ob- 
tain both, the exponent a — 11/6 as well as the corre- 
sponding prefactor. The result for u\{t) for t ^> t reads 
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ui(t) = (§ /l + 2 ?i) (*At>)~ 11/8 = 3.26856(i/ro)- 11/6 

(74) 



where we used the above results for the constants h and 
<7i . From the last equation one can see how the coupling 
between the temperature and the velocity distribution 
influences the evolution of temperature. Indeed, if there 
were no such coupling, there would be no coupling term 
in Eq. (|59|), and thus, no contribution from to the 
prefactor of Ui(t) in (|74|). This would noticeably change 
the time behavior of ui(t). On the other hand, the lead- 
ing term in the time dependence of temperature, uo(t), 
is not affected by this kind of coupling. 

On Fig. |l| and Fig. || we show the time dependence of 
the coefficient a 2 (t) of the Sonine polynomial expansion 
and of the temperature of the granular gas. The analyti- 
cal findings are compared with the numerical solution of 
the system ( ^ , ^5| ) . As it follows from the figures the ana- 
lytical theory reproduces fairly well the numerical results 
for the case of small S. 

As it follows from Fig. [I], for small 5 the following sce- 
nario of evolution of the velocity distribution takes place 
for a force-free granular gas. The initial Maxwellian dis- 
tribution evolves to a non-Maxwellian distribution, with 
the discrepancy between these two characterized by the 
second coefficient of the Sonine polynomials expansion 
02. The deviation from the Maxwellian distribution (de- 
scribed by a 2 ) quickly grows, until it saturates after a 
few collisions per particle at a "steady-state" value. At 
this instant the deviation from the Maxwellian distribu- 
tion is maximal, with the value a 2 ~ —OAS (Fig. [I], top). 
This refers to the first "fast" stage of the evolution, which 
takes place on a mean-collision time-scale ~ t c (0). Af- 
ter this maximal deviation is reached, the second "slow" 
stage of the evolution starts. At this stage a 2 decays to 
zero on "slow" time scale r ~ <5~ 1 r c (0) 3> r c (0), which 
corresponds to the time scale of the temperature evolu- 
tion (Fig. [j] middle); the decay of the coefficient a 2 (t) 
in this regime occurs according to a power law ~ t -1 / 6 
(Fig. [I], bottom). Asymptotically the Maxwellian distri- 
bution would be achieved, if the clustering process did 
not occur. 

Fig. U illustrates the significance of the first-order 
correction u\{t) in the time-evolution of temperature. 
This becomes more important as the dissipation pa- 
rameter 5 grows (Fig. || top, Fig. || middle). At 
large times the results of the first-order theory (with 
u\(t) included) practically coincide with the numer- 
ical results, while zero-order theory (without u\(t)) 
demonstrates noticeable deviations (Fig. || bottom). 
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FIG. 1. Time dependence of the second coefficient 
of the Sonine polynomial expansion azit). Time is 
given in units of mean collisional time t c (0). (Top): 
a 2 x 1000 (solid lines) for S = 0.001,0.005,0.01,0.015 
(top to bottom) together with the linear approximation 
(dashed lines); (Middle): the same as (Top) but for 
larger times; (Bottom): — <i2(t) over time (log-scale) for 
6 = 0.03, 0.01, 0.003, 0.001 (top to bottom) together with 
the power-law asymptotics ~ t~ 
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FIG. 2. Time-evolution of the reduced temperature, 
u(t) — T(t)/To. The time is given in units of mean 
collisional time t c (0). Solid line: numerical solution, 
short-dashed: ito(i) = (1 + i/ r o)~ (zero-order the- 
ory), long-dashed: u(t) — uo(t) + Sui(t) (first-order the- 
ory). (Top): for S = 0.05,0.1 (top to bottom); (Middle): 
8 — 0.15,0.25 (top to bottom); (Bottom): the same as 
(Top) but log-scale and larger ranges. 

According to our analysis of the diffusion in granular 
jas of viscoelastic particles Q , the clustering is expected 



to be retarded, compared to the case of a constant e. 
Therefore, we may assume that for the time shown on 
the figures the granular gas is still in the regime of ho- 
mogeneous cooling. 

For larger values of S the linear theory breaks down. 
Unfortunately, the equations obtained for the second 
order approximation 0(S 2 ) are too complicated to be 
treated analytically. Hence, we studied them only numer- 
ically (see Fig. ||). As compared to the case of small 5, 
an additional intermediate regime in the time-evolution 
of the velocity distribution is observed. The first "fast" 
stage of evolution takes place, as before, on the time scale 
of few collisions per particle, where maximal deviation 
from the Maxwellian distribution is achieved (Fig. |^). 
For 5 > 0.15 these maximal values of a 2 are positive. 
Then, on the second stage (intermediate regime), which 
continues 10 — 100 collisions, a 2 changes its sign and 
reaches a maximal negative deviation. Finally, on the 
third, slow stage, a 2 (t) relaxes to zero on the slow time- 
scale ~ To, just as for small 8. In Fig. [| we show the 
first stage of the time evolution of 02 (i) for systems with 
large 6. At a certain value of the dissipative parame- 
ter S the behavior changes qualitatively, i.e. the sys- 
tem then reveals another time scale as discussed above. 




FIG. 3. Time dependence of the second coeffi- 
cient of the Sonine polynomial expansion ai (t) x 100. 
Time is given in units of mean collisional time t c (0). 
5 = 0.1, 0.11, 0.12, . . . , 0.20 (bottom to top). 

Figure [| shows the numerical solution of Eqs. (44 ]45|) 
for the second Sonine coefficient a 2 (t) as a function of 
time. One can clearly distinguish the different stages of 
evolution of the velocity distribution function. 

Thus we conclude that for the case of not very small 
dissipative parameter S the time evolution of the ve- 
locity distribution function (described on the level of 
the second coefficient of the Sonine polynomials ex- 
pansion) exhibits very complicated nonmonotonic be- 
havior with few different regimes. Physically such be- 
havior is caused by existing of an additional intrin- 
sic time-scale, which describes the viscoelastic colli- 
sion and by coupling of the evolution of the velocity 
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distribution with the time-evolution of temperature. 
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FIG. 4. The second Sonine coefficient ai for S — 0.16 
over time. The numerical solutions of Eqs. (^.^) show 
all stages of evolution discussed in the text. 
The analysis which has been performed up to now has 
been addressed to the main part of the velocity dis- 
tribution function. The most important component of 
the distribution is still the Maxwellian, while deviations 



from this have been quantified in terms of the Sonine 
polynomial expansion. For very large velocities however 
this is not true and the Maxwellian distribution may not 
be used as a zero-order approximation. In the next sec- 
tion we address the problem of properties of the velocity 
distribution function for v 3> Vq. 



V. HIGH- VELOCITY TAIL OF THE 
VELOCITY-DISTRIBUTION FUNCTION 

The high- velocity tail of the velocity distribution func- 
tion in force-free granular gases was analyzed for the 
case of a constant restitution coefficient in Refs. HI). 
It was shown in these studies that for large velocities, 
c 1, the velocity distribution function behaves as 
f(c) ~ cxp(— const • c), i.e. that the tail c ^> 1 is over- 
populated, as compared to the Maxwellian distribution 
~ exp(-c 2 ). 

Here we use the scheme of analysis proposed in Rcf . [BJ . 
The same arguments as in , lead to a conclusion that 
the gain term of the collisional integral / may be ne- 
glected for c ^ 1 with respect to the loss term, which 
does not depend on the restitution coefficient. Thus, fol- 
lowing , we approximate the collision integral as 



(/,/) w-7rc/(g;t) 



(75) 



and write for c> 1 the kinetic equation ( p7| ) as 

-7rcf(c,t) (76) 
If one would use the expansion (M) (with coefficients 



a p for p > 2 discarded) to substitute it into Eq. (|76j) one 
would obtain for the second term in the left-hand side of 
© at c> 1: 



B m~3 



/Lt 2 (l + a 2 ) - \m 
5 



<Mc)S 2 (c 2 ) 



cV 



where we have used the relation 

dj/dt^a 2 4>(c)S 2 {c 2 ). 



(77) 



(78) 



a,2 according to Eq 

„2 



and the definition 



J11P of 

S 2 (c 2 ), which shows that S 2 {(?) ~ c 4 at c > 1. We 
also take into account that fj, 2 , \x± and a 2 do not depend 
on c. For the first term in the left-hand side of ([76]) and 
for the right-hand side of ( |76| ) this substitute yields cor- 
respondingly in the same limit c> 1: 



2i 

dc 
cf 



-cV 



cV 



(79) 
(80) 

From the last Eqs. @, @ and @ follows that al- 
though all terms in the Eq. (|7^) have the same factor 
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e~ c , the exponents of the power of c of the prefactor are 
different for all terms. This means inconsistency of the 
substitute (p0|), with a p for p > 2 discarded, for c 3> 1. 
Similarly, it may be shown that such inconsistency ap- 
pears for any order of the Sonine polynomial expansion. 
Indeed, using the Sonine polynomial expansion (|20| ) up 



to (arbitrary) order n, yields the estimate 
for the first term and 



„(2n+2) -c 2 



c 2n e~ 



for the second term in 
the left-hand side of (76), while for the right-hand side 
of @ one obtains - c^ 2n+1 h~ c2 . 
The exponential ansatz 



f(c,t) ~ exp{-(p(t) ■ c} 



(81) 



Note that the obtained expression ( p7\ j refers only 
for times t 3> t c (Q), when the deviations from the 
Maxwellian distribution are already well developed; it 
is not applicable for the transient times t ~ t c (0). 

As one can see from Eq. (p7|) the overpopulation (with 
respect to Maxwellian distribution) of the high-velocity 
tail decreases with time on the same time-scale ~ tq as 
a,2(t), i.e., the velocity distribution in the system ap- 
proaches the Maxwellian. However, it should be noted 
that the above considerations are valid as long as the 
overpopulation in the tail is significant to make the gain 
term in collision integral be negligible as compared to the 
loss term. 



for the kinetic equation ( |76[ ) occurs, however, to be self- 
consistent for c> 1. Substituting this into Eq. (|76| ) one 
finds that the function cp(t) in (|8l|) must satisfy 



<P + ^M2-B(/5 = nB . 



(82) 

where the time-dependence of B is given by Eq. (|35| ) 
and /i2 depends on time via a<z(t) according to (|38|). In 
linear with respect to 8 approximation ct2 ~ 8, and there- 
fore, according to (|38|), (Eo|), ^(t) ~ S'wq. Using then 
the definition ( pi| ) of 8' and expression ( |35| ) for B(t), one 
obtains 



f , 2 (t)B(t) = -r^u 3 / 5 (t) 

l2 {t) 



B(t) 



/8tt 



(83) 



with T being defined by Eq. (|6J). With Eq. (|58|) for the 
time-dependence of temperature in this approximation, 
Eq. @ reads 



TQ 



(84) 



Substituting the ansatz tp ~ (1 + t/To) u we find the ex- 
ponent, 1/ = 1/6 and the prefactor, so we arrive at the 
final result 



tp(t) = b8- L 1 + — 



1/6 



with 
6 



5 7/5 



2.25978... 



(85) 



(86) 



2 V16W 2 3 /2r(3/5) 
Thus, the velocity distribution function reads for c> 1 



f(c, t) ~ exp 



b / t 
8 V r o 



1/6' 



(87) 



VI. CONCLUSION 

We studied the velocity distribution in a homoge- 
neously cooling granular gas of viscoelastic particles 
which implies an impact velocity dependent restitution 
coefficient. We observed that contrary to the case of the 
constant restitution coefficient, the distribution function 
may not be represented in a simple scaling form, where 
the time dependence of the function occurs only via the 
time dependence of the temperature. The dependence of 
the restitution coefficient on the impact-velocity causes 
a time dependence of the coefficients of the Sonine poly- 
nomials expansion, which describes the deviation of the 
velocity distribution from the Maxwellian. 

We analyzed the time evolution of the temperature and 
of the second coefficient of the Sonine polynomials expan- 
sion ci2. Contrary to the case of the constant restitution 
coefficient, the evolution of temperature is coupled now 
with the time evolution of 0,2- 

For small values of the dissipative parameter 8 we de- 
veloped an analytical theory for the time evolution of the 
temperature of the granular gas and for the coefficient of 
the Sonine polynomials expansion the case of larger 8 
was studied numerically. We observed a complicated non- 
monotonic time-behavior of the coefficient 0,2- For small 
values of the dissipative parameter 8 we detected two dif- 
ferent stages in its time evolution: a first fast stage, which 
develops on the time scale of the mean-collision time r c , 
and the second, slow stage on the time scale ~ To ^> t c , 
on which the temperature of the granular gas changes. In 
the fast stage a maximal deviation from the Maxwellian 
distribution is achieved and then the deviation relaxes to 
zero during the second slow stage. Our numerical results 
agree well with the predictions of the analytical theory 
for small 8. 

When 8 is not small, a much more complicated time be- 
havior of the coefficient 0,2 has been revealed. In addition 
to the two stages of evolution which have been observed 
for the case of small dissipative parameter, a regime of in- 
termediate relaxation has been detected. Physically such 
complicated behavior is caused by an additional intrinsic 
time-scale, which describes the viscoelastic collision and 
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by coupling of the evolution of the velocity distribution 
with the time-evolution of temperature. 

We also analyzed the high- velocity tail of the velocity 
distribution for the case of the impact- velocity dependent 
restitution coefficient for viscoelastic particles. We found 
the same exponential overpopulation of the tail as for the 
constant restitution coefficient. However, contrary to the 
latter case where the overpopulation of the tail persists 
with time, it decreases for the impact dependent restitu- 
tion coefficient and the velocity distribution tends to the 
Maxwellian as the system evolves. 

The authors want to thank I. Goldhirsch and M. H. 
Ernst for helpful discussions. The work was supported 
by Deutsche Forschungsgemeinschaft through grant Po 
472/3-2. 



with g = V12 ■ e. Using C2 = §C 2 (Eq. (0)), one can also 



write: 



= 9 



(A5) 



We use Eq. (A5) to find the relation between the length 



of the collisional cylinders, \g\dt and \g**\dt, when the 
transformation of variables Vi*,v 2 * — > V1V2 is made in 
the collisional integral. One also needs the Jacobian of 
this transformation. To calculate this, it is convenient to 
chose the coordinate axis Z along the inter-center vector 



e, i.e., 



V\2,. 



vi, 2 



(A6) 



APPENDIX A: DERIVATION OF EQ. (17) 



Then from Eqs. (Al) follows 



The change of the particle velocities due to the inverse 
collision is described by 



V2,x] 



v i,y> 



(A7) 



tf 1 = vT-~[l + e(s*W*e 
a 2 = v 2 ** + ~[l + e(g**)]g**e 



(Al) 



where we introduce the normal relative velocity g** = 
HI2 • e and where 

e {g**) = 1 - dAa 2 ^ \g**\ 1/b + C 2 A 2 a 4/5 \g**\ 2/5 T- ■ ■ ■ 

(A2) 

according to the viscoela stic chara cter of the particles 
(see Eq. (U)). Equations (|Al[ ) and (A2) imply the con- 
servation of momentum 



and the relation: 



9 = -z{9 )9 



(A3) 
(A4) 



<z = vi, x + - {g** + g) ; 
4% = V2,z - \ {9** + g) , 

where the value of g** is expressed in te rms of g (i.e. in 
terms of v\^ z and V2, z ) according to Eq. (A5). Thus, Eqs. 
( |A7j ) explicitly express all components of the inverse- 
collision velocities v^* , v^* in terms of V\, v 2 . Straight- 
forward calculations yield for the Jacobian 



dv**dv 2 ** 



1 + ^iAa 2 / 5 |s 





1/5 



21 
25 



\g\ 2/5 T 



(A8) 



Combining (A8) with Eq. (A.5) which relates the lengths 
of collisional cylinders, one arrives at the factor Xi Eq. 
(O), in the collisional integral. 



APPENDIX II: DERIVATION OF THE MOMENTS /x P (EQS. (H©) 



To calculate the moments 

V P = -\ fdSi f dc 2 { dee (-c 12 • g) |ci 2 • el <$> (c x ) <j> (c 2 ) {l + a 2 [S 2 (c 2 ) + S 2 (c 2 )] + a\ S 2 (c 2 ) S 2 (c 2 ) } A (c? + c§) 



(1) 



it is convenient to use the center of mass velocity C and relative velocity C12 such that 



ci = C + -C12 , c 2 = C - -C12 



(2) 



The Jacobian of the transformation (0) is equal to unity and the product <j> (c\) 4> (c 2 ) transforms into 
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0(ci)0(c 2 ) 



■ exp 



3/2 



exp (-2C 2 



(2tt)3/2 — ' y 2 " 

In terms of the variables C and ci 2 the quantity [5 2 (c 2 ) + 5 2 (cf.)] in Eq. (Q) may be written as 

[s 2 (4) + s 2 (4)] = c 4 + (c ■ c 12 ) 2 + 14 + ic 2 c? 2 - 5C 2 

For 5 2 (c 2 ) 52 (c 2 ) we obtain 

5 2 (c 2 ) 5 2 (c 2 ) = if! + K 2 + K 3 + K 4 



- -4i H — 

4 12 4 



where 



7^ _ 5 6 65 4 75 2 

1 8 5 6 65 4 75 2 
K2 = 1024 Cl2 " 128 Cl2 + 128 Cl2 " 32 Cl2 

- IrV 4 j_ Ir*v 2 -t- —r 2 r 6 —r 4 r 2 —r 2 r A + ^r 2 ^ 2 

— 32 4 64 12 — g 12 32 12 16 



^4 = 7(C-ci2) -^C 2 (C-c 12 c 2 2 --C 4 (C-c 12 ) - — (C-c 12 ) c 4 2 + ^ (C-ci 2 ) 



4 1 
2 2 35 

Ci n ~~~~ 



(<5 • C12) cf 2 - y (^C? • C12 
For the quantities A (cf + cf) (p = 2, 4) we find 



1 

-( 

2 

15 



1 1 
32 



A 2 „4 , 5^ 2 



A(c 2 + cf)=-i(c 12 -e) 2 (l-e 2 ) 



and 
A 



(c 4 + c 4 ) = 2 (1 + ef (c 12 • ef (c ■ if + ± (l - £ 2 ) 2 (c 12 ■ e) 4 - i (l - , 2 ) (c 12 • ef c? 2 - & (l - c 2 ) ( 
-4(l + e) ((7-c 12 ) ((7 • e) (c 12 • g) . 
Substituting %, (§), (0) and @ into (@) and using the expansions 



(1 - e 2 ) - 2C 1 8'{t) |c 12 • el 1/5 - yC 2 <5' 2 (i) |c 12 • ef /5 + • ■ 



(1 + ef =4 



1 - CiS'(t) |c 12 • e1 1/5 + |J^' 2 (i) |c 12 • ef /5 
2/5 



4c^^(t)|c 12 - e r /o + 



one observes that fj, 2 and /i4 may be expressed in terms of the basic integrals 

Jk,i, m ,n, P , Q = J dc 12 JddJ deG (-3a ■ e) |c 12 • e\ 1+a (c 12 ) 0(tf)C*ck (c • c 12 )" 1 ((7 ■ e)" (c 12 • ej 
Namely, one has for \x 2 : 



where we define 



^0,0,0,0,2,1/5 + a 2 L ( 5) + fl 2 M (5 



— y 2 r 2 
20 d Gl 



0,0,0,2, 2/5 + a 2 L ( jrj + «2 M ( I 



— J4. 0,0, 0,2, a + ^0,0,2,0,2,0! + TF "^0,4,0, 0,2,a + 77^2,2,0, 0,2, a — 5 J2, 0,0, 0,2, a — 7 ^0,2,0,0,2, a + ~ T^0,0,0,0,2,q 

lb 2 4 4 



and 
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. , 1 T 5 65 75 1 5 

M(aj — -J8,0,0,0,2,a — - J6,0, 0,0,2, a + — ^4,0,0,0, 2, a ^"^2, 0,0, 0,2, a + ^^ 48,0,0,2,a — ^0,6, 0,0, 2, a 

65 75 3 1 1 15 

+ ^0,4, 0,0, 2, a _ 22 ^0,2, 0,0, 2, a + ^TJ ^4, 4, 0,0, 2, a + 7^6,2,0,0,2,0 + Jl, 6, 0,0, 2, a g" ^4, 2, 0,0, 2, a 

65 1 1 1 1 5 

+ Yg J%, 2, 0,0, 2, a + J Jo, 0,4, 0,2, a _ J ^2,2,2,0,2,0; _ ^J4 : 0,2 : 0,2,a ~ J 0,4,2, 0,2, a + ^ ^2, 0,2, 0,2, a 

5 T 35 T /15 X ' 



o ^0,2, 2, 0,2, a fT" J), 0,2, 0,2, a + I -jr ) Jo,0,0,0,2,a (18) 



The basic integrals may be calculated (details are given in Appendix C) and the following expressions are obtained: 

-1)p -8-2 2 
(p + a + 2)(m+l) 



w „, . '-" P - 8 ;!,~ i3p ~ [1 - (-!)«■«] r f *+%±l) r f ' + '" + ; + ° + 4 ') ("I 



for n = 0, 



(-1)p+ 1 -4-2 5 /fc + m + 4\ /f + m + p + a + 4 



Jk,l,m,l,p,a = \; +Q + 3)(m + 2) [1 - (-!)"] r ( ) r ( 77 ) (20) 



for n = 1 and 



-fc + i+p + a-l 



J *. , -^* Q = (p + a + 2)(m + l) J r 2 J F 



p+a+1 1 



m + 3 m + 1 

(21) 



for n = 2. 

When we compare Eqs. (fl6|) and (^q) for /12 and use relations (|l9|), ( p(i| ) and (|2l]) for the basic integrals, we find 
Ai = ^^0,0,0,0,2,1/5 = 2V2^2 1 / 10 r Ci = 6.48562 ... = cj (22) 



A = ^C?J ,0 1 J 0,3,2/5 = V2^2 1/6 r (y) C? - 9. 



28569 . . . = wi . (23) 



Computing from the basic integrals L{a) and M(a) for a = | and a = |, and using the relation for the Gamma- 
function T(x + 1) = xT(x), we obtain 

M = ^ (|) = |.o (24) 

*=M§)=^ (25) 

A3 = (I) = Ji^ (26) 

A = (?) = (27) 

6 20 1 V5/ 640000 v ; 

Similar calculations may be performed for /14 and yield Eq. ( |39| ) with the coefficients expressed in terms of the 
basic integrals: 

Bi = 4 (■/::„ ;i - ^0,0,0,2,2,0) = 4(V2^ - V^tt) = (28) 

^2 = —4 ( — Jo 4, 0,2, 2,0 — ^4, 0,1, 1,1,0 + <A), 0,2, 2, 2,0 — Jo, 0,3, 1,1,0 + J4.0, 0,2,2,0 — 777^0.4,1,1,1,0 + 77^2,2,0,2,2,0 

\16 16 2 

1 5 5 15 15 

— 2^2,2,1,1,1,0 — 7^0,2,0,2,2,0 + 7^0,2,1,1,1,0 — 5^2,0,0,2,2,0 + 5^2,0,1,1,1,0 + y ^0,0,0,2,2,0 7^0,0,1,1,1,0 

= 4V2^ (29) 
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^4 — 4Ci ^0,0,0,2,2,1/5 _ 2^°'°' 1 ' 1 » 1 ' 1 / 5 + Ys^ ' 2 ' ' ' 2 - 1 / 5 + ^^2,0,0,0,2,1/5^ = — \/27r2 1/10 r ^Yo) ^ = IT^ 

&7 = C\ ^ — J2, 0,0, 0,2, 2/5 + ^^0,2,0,0,2,2/5 + "g" ^0,0,0,2,2,2/5 + ^Jo, 0,0, 0,4, 2/5 ~ g ^0,0,1,1,1,2/5^ 

We do not give the expressions for the other few coefficients Bk in terms of the basic integrals, since they are too 
much cumbersome to be written explicitly. Computations of these is straightforward and yields the result: 

B 3 = gV2^ (32) 

1806 . , 

B 5 = — c (33) 

B * = liSo- ^ 

K 149054 
R 348424 

Da = Wl (36) 
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APPENDIX III: CALCULATIONS OF THE 
BASIC INTEGRALS J K ,L,M,N,p, a 

In this appendix we give some details for the calcu- 
lations of the basic integrals Jk,i,m,n,p,a- We need only 
integrals for n — 0,1,2. Evaluation of the integral for 
n = is straightforward, however, for n = 1, 2 it requires 
some tricks which are described e.g. in ]2^| . 

For n = 1 the basic integrals may be written as 



(1) 



Jk,l J dg J dC(j>(g)(j){C) x 

cV(<?-sT (<?•!(<?)) 

with <f = c*i2 and with the vectorial integral 

1(9) = J dne\g-e\ a (g ■ ef 

with the short-hand notation dfi = deO(—g ■ e) \g ■ e| 
Similarly, for n = 2 one can write 



(2) 



J dg J dC 4>(g)4>(C) x 
C k g l (C-g) m C-H(g)-C (3) 
where the dyad H(g) is given by 

H{g)= fd^eoe\g-e\ a (g-e) p , (4) 



and where o denotes direct vector product. Due to sym- 
metry one can write 1(g) = gG(g), where the function 
G(g) may be found from the equation 



g- 1(g) = g 2 G(g) = deG(-g-e) \g ■ e\ 1+a (g ■ e) 



1 + a f7t. J\P +1 



(5) 



The integral in the right-hand side of (||) may be evalu- 
ated using spherical coordinates: 



deQ(-g-e) \g ■ ef (g ■ e) r = f^^f^ • (6) 



This yields the function G(g) and, thus, the vectorial in- 
tegral 



1(g) = 2tt(-1) 



p+1 -9- 



p+a 



p + a + 3 



(7) 



For the dyad H (g) one can also use symmetry arguments 
to write 



H(g) =A(g)gog + B(g)g 2 U , 



(8) 



where U is a unit dyad (i.e. diagonal matrix). Multiply- 
ing H from both sides by g and then taking the trace, 
we obtain the set of equations for the functions A(g) and 
B( 9 ): 



g-H-g = Ag 4 + Bg 4 = 



(9) 



deQ(-g-e) \g • e\ 1+a (g ■ e) p+2 = 



p + a + 4 



.9' 



and 

TrH = Ag 2 + 3Bg 2 = 



(10) 



deO(-g-e) \g ■ e\ 1+a (g ■ gf = 

p + a + 2 
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Solving the set ©, @ for and 5(5) we obtain 
2TT(-l)PgP +a - 1 



H = 



(p + a + A){p + a + 2) 



(p + a + l)gog + g'U 



(11) 



With Eqs. (|7|) and ([11]) the basic integrals Jk,l,m,n,p,a for 
n = 1 and n = 2 can be reduced to the integrals 



dg / d<3(j>(g)(l)(C)C k+,/1 g l +P+ a +^ (C ■ g) m + v * (12) 



with v\ = 0, i>i = 0, 1S3 = 1 to evaluate the integral for 
n = 1 and with i^i = 0, 1^2 = — 1) ^3 = 2 and v\ = 2, 
^2 = 1, V3 = for n = 2. The computation of these in- 
tegrals is straightforward and yield the final result (19), 
( pp| ) and (|2l]) which has been given above. 
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